ID Analyses
2023-12-07
Calorie Data
Data Preparation
Let’s start by loading in our packages:
library(tidyverse)
library(readxl)
library(rstatix)
library(ggpubr)
library(ez)
library(emmeans)
library(car)
library(Hmisc)
library(corrplot)
library(ggsignif)
library(writexl)And let’s read our excel data in
#Vector for column names and diet names
column_names <- c("rat", "zero", "twenty", "forty", "sixty", "eighty", "hundred", "h2o", "spill", "condition", "day")
diet_names <- c("zero", "twenty", "forty", "sixty", "eighty", "hundred")
condition_names <- c("HI", "LO", "VEH", "NAL")
#Read in data
data <- read_excel("ID Test Data.xlsx", sheet = 6, range = "A3:K50", col_names = column_names)
#turn negatives to zeroes
data_filtered <- data
data_filtered[data_filtered < 0] <- 0
#filter out ID11 (because we missed on his surgery)
data_filtered <- data_filtered %>%
filter(rat != "ID11")Now we need to pivot the data into long format:
data_long <- data_filtered %>%
pivot_longer(cols = diet_names, names_to = "diet", values_to = "kcal")Also, because spillage and water consumption are in grams, for now I’m just going to drop it:
data_long_foodOnly <- data_long %>%
subset(select = -c(h2o, spill))And with that I think we have the data looking how we want it.
Basic visualizations
Now lets try visualizing things before we do any inferential analyses. For a start, lets plot average caloric intake over drug condition
#mean calories for each drug condition (collapsed across days, conditions, and rats)
meanKcals_groupedByDrug <- data_long_foodOnly %>%
group_by(condition) %>%
summarise(kcal = mean(kcal))
#SD's
SDKcals_groupedByDrug <- data_long_foodOnly %>%
group_by(condition) %>%
summarise(sd = sd(kcal),
n = n(),
se = sd/sqrt(n))
#Merge data frames
data_groupedByDrug <- inner_join(meanKcals_groupedByDrug,
SDKcals_groupedByDrug,
by = "condition" )
data_groupedByDrug$condition <- factor(data_groupedByDrug$condition,
levels = c("HI", "LO", "VEH", "NAL"))
#Make plot
my_drug_labels <- setNames(c("High DAMGO", "Low DAMGO", "Naltrexone", "Saline"),
c("HI", "LO", "NAL", "VEH"))
data_groupedByDrug %>%
ggplot(mapping = aes(x = condition,
y = kcal,
fill = condition,
ymin = kcal - se,
ymax = kcal + se)) +
geom_col(show.legend = FALSE) +
geom_errorbar(position = position_dodge(0.9),
width = 0.5,
size = 0.2) +
scale_x_discrete(labels = my_drug_labels) +
labs(x = "Drug Condition",
y = "Intake (kcal)") +
theme_classic()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
And to get a sense of between-subject variability, lets do a similar analysis grouping by rat
#mean calories for each rat (collapsed across days, conditions, and drug condition)
meanKcals_groupedByRat <- data_long_foodOnly %>%
group_by(rat) %>%
summarise(mean(kcal))
#SD's
SDKcals_groupedByRat <- data_long_foodOnly %>%
group_by(rat) %>%
summarise(sd(kcal))
#Merge data frames
data_groupedByRat <- inner_join(meanKcals_groupedByRat,
SDKcals_groupedByRat,
by = "rat" ) %>%
rename("mean_kcal" = "mean(kcal)",
"sd_kcal" = "sd(kcal)")
rat_names <- c("ID1", "ID2", "ID3", "ID4", "ID5", "ID6", "ID7", "ID8", "ID9", "ID10", "ID12")
data_groupedByRat$rat <- factor(data_groupedByRat$rat,
levels = rat_names)
#Make plot
data_groupedByRat %>%
ggplot(mapping = aes(x = rat,
y = mean_kcal,
fill = rat)) +
geom_col() +
geom_errorbar(aes(ymin = mean_kcal - sd_kcal,
ymax = mean_kcal + sd_kcal))
And just as a rough replication of the last study, lets group by diet too:
#mean calories for each rat (collapsed across days, conditions, and drug condition)
meanKcals_groupedByDiet <- data_long_foodOnly %>%
group_by(diet) %>%
summarise(kcal = mean(kcal))
#SD's
SEKcals_groupedByDiet <- data_long_foodOnly %>%
group_by(diet) %>%
summarise(sd = sd(kcal),
n = n(),
se = sd/sqrt(n))
#Merge data frames
data_groupedByDiet <- inner_join(meanKcals_groupedByDiet,
SEKcals_groupedByDiet,
by = "diet" )
data_groupedByDiet$diet <- factor(data_groupedByDiet$diet,
levels = diet_names)
#Make plot
my_diet_labels <- setNames(c("0%", "20%", "40%", "60%", "80%", "100%"),
diet_names)
data_groupedByDiet %>%
ggplot(mapping = aes(x = diet,
y = kcal,
fill = diet,
ymin = kcal - se,
ymax = kcal + se)) +
geom_col(show.legend = FALSE) +
geom_errorbar(position = position_dodge(0.9),
width = 0.5,
size = 0.2) +
scale_x_discrete(labels = my_diet_labels) +
labs(x = "Diet (% sugar by weight)",
y = "Intake (kcal)") +
theme_classic()
So here it looks like on average the eighty percent diet was a hit. However, it seems like across all of these dimensions, as overall caloric intake increases, the variation in intake also increases.
Inferential Stats
Now let’s do the hypothesis test:
res.aov = anova_test(
data = data_long_foodOnly, dv = kcal, wid = rat,
within = c(condition, diet)
)
get_anova_table(res.aov)## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 condition 3.00 30.00 29.188 4.95e-09 * 0.070
## 2 diet 2.17 21.71 7.929 2.00e-03 * 0.313
## 3 condition:diet 15.00 150.00 4.536 4.66e-07 * 0.145#trying it a different way
exp2 <- aov(kcal ~ diet*condition + Error(rat/(diet*condition)), data = data_long_foodOnly)
summary(exp2)##
## Error: rat
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 10 980.2 98.02
##
## Error: rat:diet
## Df Sum Sq Mean Sq F value Pr(>F)
## diet 5 17228 3446 7.929 1.46e-05 ***
## Residuals 50 21727 435
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: rat:condition
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 3 2855.2 951.7 29.19 4.95e-09 ***
## Residuals 30 978.2 32.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: rat:diet:condition
## Df Sum Sq Mean Sq F value Pr(>F)
## diet:condition 15 6394 426.2 4.536 4.66e-07 ***
## Residuals 150 14094 94.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#And a third way because the first two don't match
kcals.aov <- ezANOVA(data = data_long_foodOnly,
dv = .(kcal),
wid = .(rat),
within = .(diet, condition),
return_aov = TRUE,
type = 1)## Warning: Converting "rat" to factor for ANOVA.## Warning: Converting "diet" to factor for ANOVA.## Warning: Converting "condition" to factor for ANOVA.kcals.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 1 diet 5 50 7.929407 1.464019e-05 * 0.31887886
## 2 condition 3 30 29.187765 4.948465e-09 * 0.07200046
## 3 diet:condition 15 150 4.536271 4.661740e-07 * 0.14802151So these two methods of calculation slightly disagree about the p-value for the main effect of drug, but they both clearly show two main effects and an interaction. The p-values are actually suspiciously small here; we might have to circle back to this.
In any case, I want to get a more detailed visualization before I decide what post-hoc analyses I want to run:
ggboxplot(data_long_foodOnly,
x = "diet", y = "kcal",
color = "condition"
)
#mean calories for each rat (collapsed across days, conditions, and drug condition)
meanKcals_groupedByDiet <- data_long_foodOnly %>%
group_by(diet, condition) %>%
summarise(kcal = mean(kcal))## `summarise()` has grouped output by 'diet'. You can override using the
## `.groups` argument.#SD's
SEKcals_groupedByDiet <- data_long_foodOnly %>%
group_by(diet, condition) %>%
summarise(n = n(),
sd = sd(kcal),
se = sd/sqrt(n))## `summarise()` has grouped output by 'diet'. You can override using the
## `.groups` argument.#Merge data frames
data_groupedByDiet <- inner_join(meanKcals_groupedByDiet,
SEKcals_groupedByDiet,
by = c("diet", "condition"))
data_groupedByDiet$diet <- factor(data_groupedByDiet$diet,
levels = diet_names)
#Make plot
my_x_labels <- setNames(c("High DAMGO", "Low DAMGO", "Naltrexone", "Saline"),
c("HI", "LO", "NAL", "VEH"))
my_fill_labels <- setNames(c("0%", "20%", "40%", "60%", "80%", "100%"),
diet_names)
data_groupedByDiet %>%
ggplot(mapping = aes(x = condition,
y = kcal,
fill = diet,
ymin = kcal - se,
ymax = kcal + se)) +
geom_col(position = "dodge") +
geom_errorbar(position = position_dodge(0.9),
width = 0.5,
size = 0.2) +
scale_x_discrete(labels = my_x_labels) +
scale_fill_discrete(labels = my_fill_labels) +
labs(x = "Drug Condition",
y = "Intake (kcal)",
fill = "Diet (% sugar)") +
theme_classic() +
theme_classic(base_size = 16) +
theme(axis.text.x = element_text(size = 12))
So at least descriptively, this looks like good news! Consumption was highest of the sixty and eighty percent diet, and within each group, drugs look like they behaved the way we were expecting. Additionally, the significant interaction suggests that the drug effect was diet-dependent; the high dose of DAMGO spiked calorie intake a lot more for the eighty percent diet than it did for any of the others.
However, what I’m concerned about is that calculating things based on calories might actually be giving us an underestimate of what’s going on. That’s because the diets favored by the rats are actually the ones with lower caloric density – the ones with higher sugar percentages by weight. I want to check this by repeating the same analysis pipeline from the data in grams.
Weight Data
Data preparation
#Vector for column names and diet names
column_names <- c("rat", "zero", "twenty", "forty", "sixty", "eighty", "hundred", "h2o", "spill", "condition", "day")
diet_names <- c("zero", "twenty", "forty", "sixty", "eighty", "hundred")
#Read in data
data_grams <- read_excel("ID Test Data.xlsx", sheet = 5, range = "A3:K50", col_names = column_names)
#turn negatives to zeroes
data_grams_filtered <- data_grams
data_grams_filtered[data_grams_filtered < 0] <- 0
#filter out ID11
data_grams_filtered <-data_grams_filtered %>%
filter(rat != "ID11")
#Pivot long
data_grams_long <- data_grams_filtered %>%
pivot_longer(cols = diet_names, names_to = "diet", values_to = "grams")
#Drop spill and water
data_grams_long_foodOnly <- data_grams_long %>%
subset(select = -c(h2o, spill))Basic visualizations
#mean calories for each drug condition (collapsed across days, conditions, and rats)
meanGrams_groupedByDrug <- data_grams_long_foodOnly %>%
group_by(condition) %>%
summarise(grams = mean(grams))
#SD's
SDGrams_groupedByDrug <- data_grams_long_foodOnly %>%
group_by(condition) %>%
summarise(sd = sd(grams),
n = n(),
se = sd/sqrt(n))
#Merge data frames
data_grams_groupedByDrug <- inner_join(meanGrams_groupedByDrug,
SDGrams_groupedByDrug,
by = "condition" )
data_grams_groupedByDrug$condition <- factor(data_grams_groupedByDrug$condition,
levels = c("HI", "LO", "VEH", "NAL"))
#Make plot
data_grams_groupedByDrug %>%
ggplot(mapping = aes(x = condition,
y = grams,
fill = condition,
ymin = grams - se,
ymax = grams + se)) +
geom_col(show.legend = FALSE) +
geom_errorbar(position = position_dodge(0.9),
width = 0.5,
size = 0.2) +
scale_x_discrete(labels = my_drug_labels) +
labs(x = "Drug Condition",
y = "Intake (grams)") +
theme_classic()
#mean calories for each rat (collapsed across days, conditions, and drug condition)
meanGrams_groupedByRat <- data_grams_long_foodOnly %>%
group_by(rat) %>%
summarise(mean(grams))
#SD's
SDGrams_groupedByRat <- data_grams_long_foodOnly %>%
group_by(rat) %>%
summarise(sd(grams))
#Merge data frames
data_grams_groupedByRat <- inner_join(meanGrams_groupedByRat,
SDGrams_groupedByRat,
by = "rat" ) %>%
rename("mean_grams" = "mean(grams)",
"sd_grams" = "sd(grams)")
rat_names <- c("ID1", "ID2", "ID3", "ID4", "ID5", "ID6", "ID7", "ID8", "ID9", "ID10", "ID11", "ID12")
data_grams_groupedByRat$rat <- factor(data_grams_groupedByRat$rat,
levels = rat_names)
#Make plot
data_grams_groupedByRat %>%
ggplot(mapping = aes(x = rat,
y = mean_grams,
fill = rat)) +
geom_col() +
geom_errorbar(aes(ymin = mean_grams - sd_grams,
ymax = mean_grams + sd_grams))
#mean calories for each rat (collapsed across days, conditions, and drug condition)
meanGrams_groupedByDiet <- data_grams_long_foodOnly %>%
group_by(diet) %>%
summarise(grams = mean(grams))
#SD's
SDGrams_groupedByDiet <- data_grams_long_foodOnly %>%
group_by(diet) %>%
summarise(sd = sd(grams),
n = n(),
se = sd/sqrt(n))
#Merge data frames
data_grams_groupedByDiet <- inner_join(meanGrams_groupedByDiet,
SDGrams_groupedByDiet,
by = "diet" )
data_grams_groupedByDiet$diet <- factor(data_grams_groupedByDiet$diet,
levels = diet_names)
#Make plot
data_grams_groupedByDiet %>%
ggplot(mapping = aes(x = diet,
y = grams,
fill = diet,
ymin = grams - se,
ymax = grams + se)) +
geom_col(show.legend = FALSE) +
geom_errorbar(position = position_dodge(0.9),
width = 0.5,
size = 0.2) +
scale_x_discrete(labels = my_diet_labels) +
labs(x = "Diet (% sugar by weight)",
y = "Intake (grams)") +
theme_classic()
Inferential Stats
res.aovGrams = anova_test(
data = data_grams_long_foodOnly, dv = grams, wid = rat,
within = c(condition, diet)
)
get_anova_table(res.aovGrams)## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 condition 3.00 30.00 30.621 2.90e-09 * 0.080
## 2 diet 1.83 18.32 10.690 1.00e-03 * 0.367
## 3 condition:diet 15.00 150.00 5.075 4.87e-08 * 0.171#trying it a different way
exp2Grams <- aov(grams ~ diet*condition + Error(rat/(diet*condition)), data = data_grams_long_foodOnly)
summary(exp2Grams) ##
## Error: rat
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 10 28.04 2.804
##
## Error: rat:diet
## Df Sum Sq Mean Sq F value Pr(>F)
## diet 5 698.6 139.72 10.69 5.06e-07 ***
## Residuals 50 653.5 13.07
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: rat:condition
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 3 104.39 34.80 30.62 2.9e-09 ***
## Residuals 30 34.09 1.14
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: rat:diet:condition
## Df Sum Sq Mean Sq F value Pr(>F)
## diet:condition 15 247.3 16.490 5.075 4.87e-08 ***
## Residuals 150 487.4 3.249
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ggboxplot(data_grams_long_foodOnly,
x = "diet", y = "grams",
color = "condition"
)
#And a third way because the first two don't match
grams.aov <- ezANOVA(data = data_grams_long_foodOnly,
dv = .(grams),
wid = .(rat),
within = .(diet, condition),
return_aov = TRUE,
type = 1)## Warning: Converting "rat" to factor for ANOVA.## Warning: Converting "diet" to factor for ANOVA.## Warning: Converting "condition" to factor for ANOVA.grams.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 1 diet 5 50 10.690176 5.057562e-07 * 0.3728687
## 2 condition 3 30 30.621458 2.902644e-09 * 0.0815962
## 3 diet:condition 15 150 5.074863 4.872602e-08 * 0.1739020#Figure 2
#mean calories for each rat (collapsed across days, conditions, and drug condition)
meanGrams <- data_grams_long_foodOnly %>%
group_by(diet, condition) %>%
summarise(grams = mean(grams))## `summarise()` has grouped output by 'diet'. You can override using the
## `.groups` argument.#SD's
SEGrams <- data_grams_long_foodOnly %>%
group_by(diet, condition) %>%
summarise(n = n(),
sd = sd(grams),
se = sd/sqrt(n))## `summarise()` has grouped output by 'diet'. You can override using the
## `.groups` argument.#Merge data frames
data_joined <- inner_join(meanGrams,
SEGrams,
by = c("diet", "condition"))
data_joined$diet <- factor(data_joined$diet,
levels = diet_names)
#Make plot
my_x_labels <- setNames(c("High DAMGO", "Low DAMGO", "Naltrexone", "Saline"),
c("HI", "LO", "NAL", "VEH"))
my_fill_labels <- setNames(c("0%", "20%", "40%", "60%", "80%", "100%"),
diet_names)
data_joined %>%
ggplot(mapping = aes(x = condition,
y = grams,
fill = diet,
ymin = grams - se,
ymax = grams+ se)) +
geom_col(position = "dodge") +
geom_errorbar(position = position_dodge(0.9),
width = 0.5,
size = 0.2) +
scale_x_discrete(labels = my_x_labels) +
scale_fill_discrete(labels = my_fill_labels) +
labs(x = "Drug Condition",
y = "Intake (grams)",
fill = "Diet (% sugar)") +
theme_classic(base_size = 16) +
theme(axis.text.x = element_text(size = 12))
So the story in terms of grams looks just about the same as it does in terms of calories. However, there’s enough variability in the data that two possibilities are jumping into my head:
1.) Each individual rat may have very different preferences. Some may be sugar-preferrers while others are fat-preferrers 2.) Outliers may be driving some of the results we’re seeing. For instance, we had one rat, eat 23 grams of the 80% diet on one day, and another ate 17. I feel like we need to account for that somehow.
Before I try to tackle this, however, let’s get into some post-hoc analyses.
Post Hoc Analyses
The challenge throughout this section is going to be to capture what I want to see without running a million tests. To a certain extent I’m not sure that’ll be possible, but I do think I should make an effort.
Looking back at my prospectus, the questions I want to answer are:
1.) What ratio do the rats prefer (should be a replication from experiment 1)
2.) Do the rats prefer blend diets to the single nutrient ones (again, a replication)
3.) Does opioid blockade shift diet preference to lower fat or abolish preferences altogether?
4.) Does opioid stimulation shift preference to higher fat?
5.) Do opioids overall increase consumption?
6.) Does naltrexone overall decrease consumption?
Question 1
To figure out which diets the rats preferred overall, let’s first visualize the overall pattern. We actually did a similar plot earlier, but lets try a boxplot this time:
#For calories
data_long_foodOnly %>%
ggplot(mapping = aes(x = factor(diet, level = diet_names),
y = kcal, fill = diet)) +
geom_boxplot()
#For grams
data_grams_long_foodOnly %>%
ggplot(mapping = aes(x = factor(diet, level = diet_names),
y = grams, fill = diet)) +
geom_boxplot()
And now that we have a visual, we can run pairwise paired t-tests between diet conditions. To be super conservative, let’s start with Bonferroni corrections
#By calories
pairwiseKcal <- data_long_foodOnly %>%
pairwise_t_test(kcal~diet, paired = TRUE,
p.adjust.method = "bonferroni",
estimate = TRUE)
pairwiseKcal## # A tibble: 15 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 kcal eighty forty 44 44 5.66 43 1.13e-6 1.7 e-5 ****
## 2 kcal eighty hundr… 44 44 6.01 43 3.52e-7 5.28e-6 ****
## 3 kcal eighty sixty 44 44 1.82 43 7.5 e-2 1 e+0 ns
## 4 kcal eighty twenty 44 44 5.00 43 1.02e-5 1.53e-4 ***
## 5 kcal eighty zero 44 44 5.20 43 5.27e-6 7.91e-5 ****
## 6 kcal forty hundr… 44 44 1.10 43 2.76e-1 1 e+0 ns
## 7 kcal forty sixty 44 44 -3.96 43 2.79e-4 4 e-3 **
## 8 kcal forty twenty 44 44 -1.05 43 2.99e-1 1 e+0 ns
## 9 kcal forty zero 44 44 -0.986 43 3.3 e-1 1 e+0 ns
## 10 kcal hundred sixty 44 44 -4.20 43 1.31e-4 2 e-3 **
## 11 kcal hundred twenty 44 44 -1.45 43 1.55e-1 1 e+0 ns
## 12 kcal hundred zero 44 44 -1.53 43 1.34e-1 1 e+0 ns
## 13 kcal sixty twenty 44 44 3.27 43 2 e-3 3.2 e-2 *
## 14 kcal sixty zero 44 44 3.15 43 3 e-3 4.4 e-2 *
## 15 kcal twenty zero 44 44 0.0510 43 9.6 e-1 1 e+0 ns#By grams
pairwiseGrams <- data_grams_long_foodOnly %>%
pairwise_t_test(grams~diet, paired = TRUE,
p.adjust.method = "bonferroni")
pairwiseGrams## # A tibble: 15 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 grams eighty forty 44 44 5.90 43 5.16e-7 7.74e-6 ****
## 2 grams eighty hundr… 44 44 5.88 43 5.46e-7 8.19e-6 ****
## 3 grams eighty sixty 44 44 2.42 43 2 e-2 2.96e-1 ns
## 4 grams eighty twenty 44 44 5.57 43 1.52e-6 2.28e-5 ****
## 5 grams eighty zero 44 44 5.80 43 7.25e-7 1.09e-5 ****
## 6 grams forty hundr… 44 44 -0.0743 43 9.41e-1 1 e+0 ns
## 7 grams forty sixty 44 44 -4.08 43 1.9 e-4 3 e-3 **
## 8 grams forty twenty 44 44 -0.754 43 4.55e-1 1 e+0 ns
## 9 grams forty zero 44 44 -0.357 43 7.23e-1 1 e+0 ns
## 10 grams hundred sixty 44 44 -3.88 43 3.57e-4 5 e-3 **
## 11 grams hundred twenty 44 44 -0.401 43 6.9 e-1 1 e+0 ns
## 12 grams hundred zero 44 44 -0.156 43 8.77e-1 1 e+0 ns
## 13 grams sixty twenty 44 44 3.67 43 6.62e-4 1 e-2 **
## 14 grams sixty zero 44 44 3.68 43 6.4 e-4 1 e-2 **
## 15 grams twenty zero 44 44 0.283 43 7.79e-1 1 e+0 nsSo even being as conservative as possible (which we may decide is what we want) the 80% diet is significantly more consumed than all other diets except the 60%, and the 60% diet is significantly preferred over the 0%, 40%, and 100% diets.There apparently isn’t a significant difference between the 20% and 60% diets, however, which actually makes me want to take a look at the effect sizes we’re dealing with:
#By calories
data_long_foodOnly %>%
cohens_d(kcal~diet, paired = TRUE)## # A tibble: 15 × 7
## .y. group1 group2 effsize n1 n2 magnitude
## * <chr> <chr> <chr> <dbl> <int> <int> <ord>
## 1 kcal eighty forty 0.854 44 44 large
## 2 kcal eighty hundred 0.906 44 44 large
## 3 kcal eighty sixty 0.275 44 44 small
## 4 kcal eighty twenty 0.753 44 44 moderate
## 5 kcal eighty zero 0.784 44 44 moderate
## 6 kcal forty hundred 0.166 44 44 negligible
## 7 kcal forty sixty -0.597 44 44 moderate
## 8 kcal forty twenty -0.158 44 44 negligible
## 9 kcal forty zero -0.149 44 44 negligible
## 10 kcal hundred sixty -0.633 44 44 moderate
## 11 kcal hundred twenty -0.218 44 44 small
## 12 kcal hundred zero -0.230 44 44 small
## 13 kcal sixty twenty 0.493 44 44 small
## 14 kcal sixty zero 0.475 44 44 small
## 15 kcal twenty zero 0.00768 44 44 negligible#By grams
data_grams_long_foodOnly %>%
cohens_d(grams~diet, paired = TRUE)## # A tibble: 15 × 7
## .y. group1 group2 effsize n1 n2 magnitude
## * <chr> <chr> <chr> <dbl> <int> <int> <ord>
## 1 grams eighty forty 0.889 44 44 large
## 2 grams eighty hundred 0.886 44 44 large
## 3 grams eighty sixty 0.365 44 44 small
## 4 grams eighty twenty 0.840 44 44 large
## 5 grams eighty zero 0.874 44 44 large
## 6 grams forty hundred -0.0112 44 44 negligible
## 7 grams forty sixty -0.615 44 44 moderate
## 8 grams forty twenty -0.114 44 44 negligible
## 9 grams forty zero -0.0538 44 44 negligible
## 10 grams hundred sixty -0.584 44 44 moderate
## 11 grams hundred twenty -0.0605 44 44 negligible
## 12 grams hundred zero -0.0236 44 44 negligible
## 13 grams sixty twenty 0.553 44 44 moderate
## 14 grams sixty zero 0.555 44 44 moderate
## 15 grams twenty zero 0.0426 44 44 negligibleSo they’re medium to large, which is encouraging to see. Based on these analyses, I feel pretty comfortable concluding that the 80% diet was the favorite on average.
Question 2
Based on the plots we’ve generated, it’s pretty clear that the rats preferred blend diets on average over single-nutrient ones. However, just to be above reproach I want to do a formal comparison. Time for fun with contrasts! Assigning a coefficient of 1/4 for all the blend diets and a coefficient of -1/2 for the single diet conditions should give us what we want.
First I have to do some R trickery to force the software to handle my variables the way I want it to:
#kcals
data_long_foodOnly_factor <- data_long_foodOnly
names <- c(1:4)
data_long_foodOnly_factor[,names] <- lapply(data_long_foodOnly_factor[,names], factor)
data_long_foodOnly_factor$rat <- factor(data_long_foodOnly_factor$rat, levels = rat_names)
data_long_foodOnly_factor$diet <- factor(data_long_foodOnly_factor$diet, levels = diet_names)
data_long_foodOnly_factor$condition <- factor(data_long_foodOnly_factor$condition, levels = condition_names)
str(data_long_foodOnly_factor)## tibble [264 × 5] (S3: tbl_df/tbl/data.frame)
## $ rat : Factor w/ 12 levels "ID1","ID2","ID3",..: 1 1 1 1 1 1 2 2 2 2 ...
## $ condition: Factor w/ 4 levels "HI","LO","VEH",..: 4 4 4 4 4 4 2 2 2 2 ...
## $ day : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ diet : Factor w/ 6 levels "zero","twenty",..: 1 2 3 4 5 6 1 2 3 4 ...
## $ kcal : num [1:264] 0.45 1.76 3.92 6.3 10.35 ...#grams
data_grams_long_foodOnly_factor <- data_grams_long_foodOnly
names <- c(1:4)
data_grams_long_foodOnly_factor[,names] <- lapply(data_grams_long_foodOnly_factor[,names], factor)
data_grams_long_foodOnly_factor$rat <- factor(data_grams_long_foodOnly_factor$rat, levels = rat_names)
data_grams_long_foodOnly_factor$diet <- factor(data_grams_long_foodOnly_factor$diet, levels = diet_names)
data_grams_long_foodOnly_factor$condition <- factor(data_grams_long_foodOnly_factor$condition, levels = condition_names)
str(data_grams_long_foodOnly_factor)## tibble [264 × 5] (S3: tbl_df/tbl/data.frame)
## $ rat : Factor w/ 12 levels "ID1","ID2","ID3",..: 1 1 1 1 1 1 2 2 2 2 ...
## $ condition: Factor w/ 4 levels "HI","LO","VEH",..: 4 4 4 4 4 4 2 2 2 2 ...
## $ day : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ diet : Factor w/ 6 levels "zero","twenty",..: 1 2 3 4 5 6 1 2 3 4 ...
## $ grams : num [1:264] 0.05 0.22 0.56 1.05 2.07 ...And now I’m going to give it a shot:
#Set contrast coefficient matrix
#By kcals
options(contrasts = c("contr.sum", "contr.poly"))
kcal.aov <- ezANOVA(data = data_long_foodOnly_factor,
dv = .(kcal),
wid = .(rat),
within = .(diet, condition),
return_aov = TRUE,
type = 1)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.kcal.em <- emmeans(kcal.aov$aov, ~diet * condition)
contrast_names <- list(c(2, -1, -1, -1, -1, 2,
2, -1, -1, -1, -1, 2,
2, -1, -1, -1, -1, 2,
2, -1, -1, -1, -1, 2))
contrast(kcal.em,
method = contrast_names)## contrast estimate SE df
## c(2, -1, -1, -1, -1, 2, 2, -1, -1, -1, -1, 2, 2, -1, -1, -1, -134 43.5 50
## t.ratio p.value
## -3.083 0.0033#By grams
options(contrasts = c("contr.sum", "contr.poly"))
grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor,
dv = .(grams),
wid = .(rat),
within = .(diet, condition),
return_aov = TRUE,
type = 1)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.grams.em <- emmeans(grams.aov$aov, ~diet * condition)
contrast_names <- list(c(2, -1, -1, -1, -1, 2,
2, -1, -1, -1, -1, 2,
2, -1, -1, -1, -1, 2,
2, -1, -1, -1, -1, 2))
contrast(grams.em,
method = contrast_names)## contrast estimate SE df
## c(2, -1, -1, -1, -1, 2, 2, -1, -1, -1, -1, 2, 2, -1, -1, -1, -25.3 7.55 50
## t.ratio p.value
## -3.348 0.0016So if I did this right, then we did find the effect we were looking for. However, I’m a little unsure of this, so I may double-check in SPSS just for peace of mind.
Question 3
To go about answering this question, we can start by filtering the data to only include the naltrexone and saline groups:
#kcals
data_long_foodOnly_factor_q3 <- data_long_foodOnly_factor %>%
filter(condition == "VEH" | condition == "NAL")
#grams
data_grams_long_foodOnly_factor_q3 <- data_grams_long_foodOnly_factor %>%
filter(condition == "VEH" | condition == "NAL")And now we can visualize:
#kcals
ggboxplot(data_long_foodOnly_factor_q3,
x = "condition", y = "kcal",
color = "diet"
)
#grams
ggboxplot(data_grams_long_foodOnly_factor_q3,
x = "condition", y = "grams",
color = "diet"
)
So just visually, it looks like the trend of consumption rising as sugar content increases holds for both the vehicle and naltrexone groups.
Since we have a significant omnibus effect, I feel comfortable proceeding with looking at significant one-way effects within each drug category. We can start by filtering to create datasets for each specific drug group:
## kcals
# VEH dataset
data_long_foodOnly_factor_veh <- data_long_foodOnly_factor %>%
filter(condition == "VEH")
#NAL dataset
data_long_foodOnly_factor_nal <- data_long_foodOnly_factor %>%
filter(condition == "NAL")
##grams
# VEH dataset
data_grams_long_foodOnly_factor_veh <- data_grams_long_foodOnly_factor %>%
filter(condition == "VEH")
#NAL dataset
data_grams_long_foodOnly_factor_nal <- data_grams_long_foodOnly_factor %>%
filter(condition == "NAL")Now we can run one-way ANOVA’s on these data:
## kcals
#VEH
veh.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_veh,
dv = .(kcal),
wid = .(rat),
within = .(diet),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.veh.kcal.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 diet 5 50 8.839876 4.603764e-06 * 0.434701#NAL
nal.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_nal,
dv = .(kcal),
wid = .(rat),
within = .(diet),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.nal.kcal.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 diet 5 50 2.328814 0.05607767 0.1798096## grams
veh.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_veh,
dv = .(grams),
wid = .(rat),
within = .(diet),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.veh.grams.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 diet 5 50 12.29994 8.511631e-08 * 0.5238583#NAL
nal.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_nal,
dv = .(grams),
wid = .(rat),
within = .(diet),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.nal.grams.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 diet 5 50 3.677905 0.006531315 * 0.2582462So we have an omnibus effect of diet on consumption for the vehicle condition, but not the naltrexone group. Within the vehicle condition, we can proceed with pairwise comparisons:
###VEH
#Significance
data_long_foodOnly_factor_veh %>%
pairwise_t_test(kcal~diet, paired = TRUE,
p.adjust.method = "bonferroni",
estimate = TRUE)## # A tibble: 15 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 kcal zero twenty 11 11 -0.695 10 0.503 1 ns
## 2 kcal zero forty 11 11 -0.704 10 0.498 1 ns
## 3 kcal zero sixty 11 11 -2.10 10 0.062 0.928 ns
## 4 kcal zero eighty 11 11 -4.29 10 0.002 0.024 *
## 5 kcal zero hundred 11 11 2.06 10 0.066 0.996 ns
## 6 kcal twenty forty 11 11 0.470 10 0.648 1 ns
## 7 kcal twenty sixty 11 11 -1.21 10 0.253 1 ns
## 8 kcal twenty eighty 11 11 -3.97 10 0.003 0.04 *
## 9 kcal twenty hundred 11 11 1.34 10 0.21 1 ns
## 10 kcal forty sixty 11 11 -1.90 10 0.087 1 ns
## 11 kcal forty eighty 11 11 -4.07 10 0.002 0.034 *
## 12 kcal forty hundred 11 11 2.64 10 0.025 0.37 ns
## 13 kcal sixty eighty 11 11 -1.69 10 0.122 1 ns
## 14 kcal sixty hundred 11 11 2.78 10 0.019 0.29 ns
## 15 kcal eighty hundred 11 11 5.26 10 0.000366 0.005 **#Effect size
data_long_foodOnly_factor_veh %>%
cohens_d(kcal~diet, paired = TRUE)## # A tibble: 15 × 7
## .y. group1 group2 effsize n1 n2 magnitude
## * <chr> <chr> <chr> <dbl> <int> <int> <ord>
## 1 kcal zero twenty -0.209 11 11 small
## 2 kcal zero forty -0.212 11 11 small
## 3 kcal zero sixty -0.634 11 11 moderate
## 4 kcal zero eighty -1.29 11 11 large
## 5 kcal zero hundred 0.621 11 11 moderate
## 6 kcal twenty forty 0.142 11 11 negligible
## 7 kcal twenty sixty -0.366 11 11 small
## 8 kcal twenty eighty -1.20 11 11 large
## 9 kcal twenty hundred 0.404 11 11 small
## 10 kcal forty sixty -0.572 11 11 moderate
## 11 kcal forty eighty -1.23 11 11 large
## 12 kcal forty hundred 0.796 11 11 moderate
## 13 kcal sixty eighty -0.509 11 11 moderate
## 14 kcal sixty hundred 0.839 11 11 large
## 15 kcal eighty hundred 1.59 11 11 large###NAL
##kcals
#Significance
data_long_foodOnly_factor_nal %>%
pairwise_t_test(kcal~diet, paired = TRUE,
p.adjust.method = "bonferroni",
estimate = TRUE)## # A tibble: 15 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 kcal zero twenty 11 11 0.249 10 0.808 1 ns
## 2 kcal zero forty 11 11 1.13 10 0.287 1 ns
## 3 kcal zero sixty 11 11 -0.782 10 0.453 1 ns
## 4 kcal zero eighty 11 11 -1.59 10 0.143 1 ns
## 5 kcal zero hundred 11 11 1.46 10 0.174 1 ns
## 6 kcal twenty forty 11 11 1.01 10 0.337 1 ns
## 7 kcal twenty sixty 11 11 -0.977 10 0.352 1 ns
## 8 kcal twenty eighty 11 11 -2.04 10 0.069 1 ns
## 9 kcal twenty hundred 11 11 1.21 10 0.253 1 ns
## 10 kcal forty sixty 11 11 -1.60 10 0.142 1 ns
## 11 kcal forty eighty 11 11 -3.89 10 0.003 0.045 *
## 12 kcal forty hundred 11 11 1.58 10 0.146 1 ns
## 13 kcal sixty eighty 11 11 -0.340 10 0.741 1 ns
## 14 kcal sixty hundred 11 11 1.81 10 0.101 1 ns
## 15 kcal eighty hundred 11 11 4.47 10 0.001 0.018 *#Effect size
data_long_foodOnly_factor_nal %>%
cohens_d(kcal~diet, paired = TRUE)## # A tibble: 15 × 7
## .y. group1 group2 effsize n1 n2 magnitude
## * <chr> <chr> <chr> <dbl> <int> <int> <ord>
## 1 kcal zero twenty 0.0751 11 11 negligible
## 2 kcal zero forty 0.339 11 11 small
## 3 kcal zero sixty -0.236 11 11 small
## 4 kcal zero eighty -0.479 11 11 small
## 5 kcal zero hundred 0.441 11 11 small
## 6 kcal twenty forty 0.304 11 11 small
## 7 kcal twenty sixty -0.295 11 11 small
## 8 kcal twenty eighty -0.614 11 11 moderate
## 9 kcal twenty hundred 0.366 11 11 small
## 10 kcal forty sixty -0.481 11 11 small
## 11 kcal forty eighty -1.17 11 11 large
## 12 kcal forty hundred 0.475 11 11 small
## 13 kcal sixty eighty -0.103 11 11 negligible
## 14 kcal sixty hundred 0.545 11 11 moderate
## 15 kcal eighty hundred 1.35 11 11 large##grams
#Significance
data_grams_long_foodOnly_factor_nal %>%
pairwise_t_test(grams~diet, paired = TRUE,
p.adjust.method = "bonferroni",
estimate = TRUE)## # A tibble: 15 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 grams zero twenty 11 11 0.125 10 0.903 1 ns
## 2 grams zero forty 11 11 0.994 10 0.344 1 ns
## 3 grams zero sixty 11 11 -1.12 10 0.29 1 ns
## 4 grams zero eighty 11 11 -2.64 10 0.025 0.374 ns
## 5 grams zero hundred 11 11 1.33 10 0.213 1 ns
## 6 grams twenty forty 11 11 0.947 10 0.366 1 ns
## 7 grams twenty sixty 11 11 -1.19 10 0.26 1 ns
## 8 grams twenty eighty 11 11 -2.89 10 0.016 0.243 ns
## 9 grams twenty hundred 11 11 1.07 10 0.31 1 ns
## 10 grams forty sixty 11 11 -1.64 10 0.132 1 ns
## 11 grams forty eighty 11 11 -4.11 10 0.002 0.032 *
## 12 grams forty hundred 11 11 1.03 10 0.329 1 ns
## 13 grams sixty eighty 11 11 -0.666 10 0.521 1 ns
## 14 grams sixty hundred 11 11 1.76 10 0.109 1 ns
## 15 grams eighty hundred 11 11 4.42 10 0.001 0.019 *So overall, it appears that the sugarier diets are favored in the vehicle group, but these preferences are abolished (or at least are drastically reduced) in the naltrexone condition.
Question 4
I think I can repeat this pipeline to compare what happens between the DAMGO groups and the vehicle group.
Always a good idea to start with visualization:
# kcals
data_long_foodOnly_factor_q3 <- data_long_foodOnly_factor %>%
filter(condition == "HI" | condition == "LO" | condition == "VEH")
ggboxplot(data_long_foodOnly_factor_q3,
x = "condition", y = "kcal",
color = "diet"
)
# grams
data_grams_long_foodOnly_factor_q3 <- data_grams_long_foodOnly_factor %>%
filter(condition == "HI" | condition == "LO" | condition == "VEH")
ggboxplot(data_grams_long_foodOnly_factor_q3,
x = "condition", y = "grams",
color = "diet"
)
Again, the significant omnibus effect allows us to run simple one-way ANOVAs for the HI and LO dose groups:
###Filter out HI and LO groups
#kcals
# HI dataset
data_long_foodOnly_factor_hi <- data_long_foodOnly_factor %>%
filter(condition == "HI")
#LO dataset
data_long_foodOnly_factor_lo <- data_long_foodOnly_factor %>%
filter(condition == "LO")
##ANOVAs
#HI
hi.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_hi,
dv = .(kcal),
wid = .(rat),
within = .(diet),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.hi.kcal.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 diet 5 50 8.643261 5.8861e-06 * 0.4546631#LO
lo.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_lo,
dv = .(kcal),
wid = .(rat),
within = .(diet),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.lo.kcal.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 diet 5 50 4.108456 0.003339617 * 0.2793328#VEH (just so I don't have to scroll back up)
veh.kcal.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 diet 5 50 8.839876 4.603764e-06 * 0.434701## grams
# HI dataset
data_grams_long_foodOnly_factor_hi <- data_grams_long_foodOnly_factor %>%
filter(condition == "HI")
#LO dataset
data_grams_long_foodOnly_factor_lo <- data_grams_long_foodOnly_factor %>%
filter(condition == "LO")
##ANOVAs
#HI
hi.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_hi,
dv = .(grams),
wid = .(rat),
within = .(diet),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.hi.grams.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 diet 5 50 10.34611 7.522938e-07 * 0.4985924#LO
lo.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_lo,
dv = .(grams),
wid = .(rat),
within = .(diet),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.lo.grams.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 diet 5 50 4.845566 0.001087098 * 0.3126305#VEH (just so I don't have to scroll back up)
veh.grams.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 diet 5 50 12.29994 8.511631e-08 * 0.5238583So we get significant omnibus diet effects for all 3 drug groups (even though the omnibus effect for the low dose is just barely significant). Let’s run pairwise comparisons now:
## HI
#Significance
data_long_foodOnly_factor_hi %>%
pairwise_t_test(kcal~diet, paired = TRUE,
p.adjust.method = "bonferroni",
estimate = TRUE)## # A tibble: 15 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 kcal zero twenty 11 11 -0.454 10 0.66 1 ns
## 2 kcal zero forty 11 11 -0.362 10 0.725 1 ns
## 3 kcal zero sixty 11 11 -2.10 10 0.062 0.924 ns
## 4 kcal zero eighty 11 11 -4.18 10 0.002 0.028 *
## 5 kcal zero hundred 11 11 -0.0590 10 0.954 1 ns
## 6 kcal twenty forty 11 11 0.409 10 0.691 1 ns
## 7 kcal twenty sixty 11 11 -1.90 10 0.087 1 ns
## 8 kcal twenty eighty 11 11 -3.30 10 0.008 0.119 ns
## 9 kcal twenty hundred 11 11 0.452 10 0.661 1 ns
## 10 kcal forty sixty 11 11 -2.31 10 0.044 0.658 ns
## 11 kcal forty eighty 11 11 -3.83 10 0.003 0.05 *
## 12 kcal forty hundred 11 11 0.350 10 0.733 1 ns
## 13 kcal sixty eighty 11 11 -1.73 10 0.115 1 ns
## 14 kcal sixty hundred 11 11 2.14 10 0.058 0.876 ns
## 15 kcal eighty hundred 11 11 3.95 10 0.003 0.041 *#Effect size
data_long_foodOnly_factor_hi %>%
cohens_d(kcal~diet, paired = TRUE)## # A tibble: 15 × 7
## .y. group1 group2 effsize n1 n2 magnitude
## * <chr> <chr> <chr> <dbl> <int> <int> <ord>
## 1 kcal zero twenty -0.137 11 11 negligible
## 2 kcal zero forty -0.109 11 11 negligible
## 3 kcal zero sixty -0.635 11 11 moderate
## 4 kcal zero eighty -1.26 11 11 large
## 5 kcal zero hundred -0.0178 11 11 negligible
## 6 kcal twenty forty 0.123 11 11 negligible
## 7 kcal twenty sixty -0.573 11 11 moderate
## 8 kcal twenty eighty -0.996 11 11 large
## 9 kcal twenty hundred 0.136 11 11 negligible
## 10 kcal forty sixty -0.695 11 11 moderate
## 11 kcal forty eighty -1.16 11 11 large
## 12 kcal forty hundred 0.106 11 11 negligible
## 13 kcal sixty eighty -0.520 11 11 moderate
## 14 kcal sixty hundred 0.644 11 11 moderate
## 15 kcal eighty hundred 1.19 11 11 large## LO
#Significance
data_long_foodOnly_factor_lo %>%
pairwise_t_test(kcal~diet, paired = TRUE,
p.adjust.method = "bonferroni",
estimate = TRUE)## # A tibble: 15 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 kcal zero twenty 11 11 0.737 10 0.478 1 ns
## 2 kcal zero forty 11 11 1.61 10 0.139 1 ns
## 3 kcal zero sixty 11 11 -1.70 10 0.12 1 ns
## 4 kcal zero eighty 11 11 -2.24 10 0.049 0.734 ns
## 5 kcal zero hundred 11 11 0.684 10 0.51 1 ns
## 6 kcal twenty forty 11 11 0.187 10 0.855 1 ns
## 7 kcal twenty sixty 11 11 -2.21 10 0.051 0.771 ns
## 8 kcal twenty eighty 11 11 -3.01 10 0.013 0.198 ns
## 9 kcal twenty hundred 11 11 0.0329 10 0.974 1 ns
## 10 kcal forty sixty 11 11 -2.28 10 0.046 0.693 ns
## 11 kcal forty eighty 11 11 -2.90 10 0.016 0.236 ns
## 12 kcal forty hundred 11 11 -0.141 10 0.891 1 ns
## 13 kcal sixty eighty 11 11 -0.0354 10 0.972 1 ns
## 14 kcal sixty hundred 11 11 2.30 10 0.044 0.661 ns
## 15 kcal eighty hundred 11 11 2.92 10 0.015 0.229 ns#Effect size
data_long_foodOnly_factor_lo %>%
cohens_d(kcal~diet, paired = TRUE)## # A tibble: 15 × 7
## .y. group1 group2 effsize n1 n2 magnitude
## * <chr> <chr> <chr> <dbl> <int> <int> <ord>
## 1 kcal zero twenty 0.222 11 11 small
## 2 kcal zero forty 0.484 11 11 small
## 3 kcal zero sixty -0.513 11 11 moderate
## 4 kcal zero eighty -0.676 11 11 moderate
## 5 kcal zero hundred 0.206 11 11 small
## 6 kcal twenty forty 0.0564 11 11 negligible
## 7 kcal twenty sixty -0.667 11 11 moderate
## 8 kcal twenty eighty -0.906 11 11 large
## 9 kcal twenty hundred 0.00990 11 11 negligible
## 10 kcal forty sixty -0.686 11 11 moderate
## 11 kcal forty eighty -0.875 11 11 large
## 12 kcal forty hundred -0.0426 11 11 negligible
## 13 kcal sixty eighty -0.0107 11 11 negligible
## 14 kcal sixty hundred 0.694 11 11 moderate
## 15 kcal eighty hundred 0.881 11 11 large## VEH
#Significance
data_long_foodOnly_factor_veh %>%
pairwise_t_test(kcal~diet, paired = TRUE,
p.adjust.method = "bonferroni",
estimate = TRUE)## # A tibble: 15 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 kcal zero twenty 11 11 -0.695 10 0.503 1 ns
## 2 kcal zero forty 11 11 -0.704 10 0.498 1 ns
## 3 kcal zero sixty 11 11 -2.10 10 0.062 0.928 ns
## 4 kcal zero eighty 11 11 -4.29 10 0.002 0.024 *
## 5 kcal zero hundred 11 11 2.06 10 0.066 0.996 ns
## 6 kcal twenty forty 11 11 0.470 10 0.648 1 ns
## 7 kcal twenty sixty 11 11 -1.21 10 0.253 1 ns
## 8 kcal twenty eighty 11 11 -3.97 10 0.003 0.04 *
## 9 kcal twenty hundred 11 11 1.34 10 0.21 1 ns
## 10 kcal forty sixty 11 11 -1.90 10 0.087 1 ns
## 11 kcal forty eighty 11 11 -4.07 10 0.002 0.034 *
## 12 kcal forty hundred 11 11 2.64 10 0.025 0.37 ns
## 13 kcal sixty eighty 11 11 -1.69 10 0.122 1 ns
## 14 kcal sixty hundred 11 11 2.78 10 0.019 0.29 ns
## 15 kcal eighty hundred 11 11 5.26 10 0.000366 0.005 **#Effect size
data_long_foodOnly_factor_veh %>%
cohens_d(kcal~diet, paired = TRUE)## # A tibble: 15 × 7
## .y. group1 group2 effsize n1 n2 magnitude
## * <chr> <chr> <chr> <dbl> <int> <int> <ord>
## 1 kcal zero twenty -0.209 11 11 small
## 2 kcal zero forty -0.212 11 11 small
## 3 kcal zero sixty -0.634 11 11 moderate
## 4 kcal zero eighty -1.29 11 11 large
## 5 kcal zero hundred 0.621 11 11 moderate
## 6 kcal twenty forty 0.142 11 11 negligible
## 7 kcal twenty sixty -0.366 11 11 small
## 8 kcal twenty eighty -1.20 11 11 large
## 9 kcal twenty hundred 0.404 11 11 small
## 10 kcal forty sixty -0.572 11 11 moderate
## 11 kcal forty eighty -1.23 11 11 large
## 12 kcal forty hundred 0.796 11 11 moderate
## 13 kcal sixty eighty -0.509 11 11 moderate
## 14 kcal sixty hundred 0.839 11 11 large
## 15 kcal eighty hundred 1.59 11 11 largeSo without this getting super tedious, it looks like eighty percent takes it across all the two DAMGO conditions as well.
Questions 5 & 6
Keeping in mind that we have a significant main effect of drug on consumption, we can graph consumption of each diet:
data_long_foodOnly %>%
ggplot(mapping = aes(x = factor(condition, level = condition_names),
y = kcal, fill = condition)) +
geom_boxplot()
And we can then do pairwise paired t-tests:
## Kcal
# Significance
data_long_foodOnly %>%
pairwise_t_test(kcal~condition, paired = TRUE,
p.adjust.method = "bonferroni",
estimate = TRUE)## # A tibble: 6 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 kcal HI LO 66 66 1.40 65 0.167 1 ns
## 2 kcal HI NAL 66 66 2.93 65 0.005 0.028 *
## 3 kcal HI VEH 66 66 2.67 65 0.01 0.058 ns
## 4 kcal LO NAL 66 66 2.04 65 0.046 0.274 ns
## 5 kcal LO VEH 66 66 1.82 65 0.073 0.439 ns
## 6 kcal NAL VEH 66 66 -0.760 65 0.45 1 ns#Effect size
data_long_foodOnly %>%
cohens_d(kcal~condition, paired = TRUE)## # A tibble: 6 × 7
## .y. group1 group2 effsize n1 n2 magnitude
## * <chr> <chr> <chr> <dbl> <int> <int> <ord>
## 1 kcal HI LO 0.172 66 66 negligible
## 2 kcal HI NAL 0.360 66 66 small
## 3 kcal HI VEH 0.328 66 66 small
## 4 kcal LO NAL 0.251 66 66 small
## 5 kcal LO VEH 0.224 66 66 small
## 6 kcal NAL VEH -0.0935 66 66 negligible## Grams
data_grams_long_foodOnly %>%
pairwise_t_test(grams~condition, paired = TRUE,
p.adjust.method = "bonferroni",
estimate = TRUE)## # A tibble: 6 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 grams HI LO 66 66 1.52 65 0.132 0.792 ns
## 2 grams HI NAL 66 66 3.03 65 0.004 0.021 *
## 3 grams HI VEH 66 66 2.69 65 0.009 0.054 ns
## 4 grams LO NAL 66 66 2.19 65 0.032 0.192 ns
## 5 grams LO VEH 66 66 1.88 65 0.064 0.386 ns
## 6 grams NAL VEH 66 66 -1.00 65 0.319 1 ns#Effect size
data_grams_long_foodOnly %>%
cohens_d(grams~condition, paired = TRUE)## # A tibble: 6 × 7
## .y. group1 group2 effsize n1 n2 magnitude
## * <chr> <chr> <chr> <dbl> <int> <int> <ord>
## 1 grams HI LO 0.188 66 66 negligible
## 2 grams HI NAL 0.372 66 66 small
## 3 grams HI VEH 0.332 66 66 small
## 4 grams LO NAL 0.270 66 66 small
## 5 grams LO VEH 0.232 66 66 small
## 6 grams NAL VEH -0.124 66 66 negligibleSo it looks like the high dose did increase consumption over the vehicle dose, but the low dose didn’t, and naltrexone didn’t significantly lower it. If anything, it looks like drug didn’t lead to huge group level changes, but more just affected the amount of variability across groups.
Additional analysis
We’ve seen so far that drug mildly affects the average amount consumed in general, and that DAMGO exacerbates the preference for the 80% diet while naltrexone flattens it. However, one possibility I want to explore is whether different things are happening for animals that prefer fat vs. those that prefer sugar. This suggests two questions:
7.) Are there differences in baseline preference between animals, and if so, how do we determine those?
8.) Does the pattern of data we obtain differ between fat and sugar preferrers?
In my prospectus, I mention a paper (Gosnell et al., 1990), that describes a procedure for relating baseline preferences to drug effects. I’m going to try to follow their procedure here.
Question 7
As a super simple first pass, I think I can start by filtering the calorie data such that we’re only looking at the vehicle condition and we’re only including the two pure macronutrient conditions:
data_long_foodOnly_pure <- data_long_foodOnly %>%
filter(condition == "VEH") %>%
filter(diet == "zero" | diet == "hundred")
data_long_foodOnly_pure %>%
ggplot(mapping = aes(x = rat,
y = kcal,
fill = diet)) +
geom_col(position = "dodge")
And it’s only at this point that I realize the kcal data probably isn’t the best to be looking at here, because fat is more the twice as caloric as sugar. Let’s repeat the analysis with weight data:
data_grams_long_foodOnly_pure <- data_grams_long_foodOnly %>%
filter(condition == "VEH") %>%
filter(diet == "zero" | diet == "hundred")
data_grams_long_foodOnly_pure %>%
ggplot(mapping = aes(x = rat,
y = grams,
fill = diet)) +
geom_col(position = "dodge")
This still feels like it’s giving us an incomplete picture though. It actually looks like the Gosnell paper had dedicated days from which they detemined baseline preferences. For our data, I can see a few possibilities for which days to use:
1.) Some average of the diet acclimation period. The benefit is that the rats are totally injection naive at this point. Con is that since rats weren’t acclimated to diets in the same order, relative experience in the chambers between fat and sugar diets varies across rats.
2.) An average of the training days where rats were given all 6 diets
3.) The maintenance days in between injections
My gut is telling me that the training days with all 6 diets is the best way to go, so I’m going to start down that road and see what happens.
Let’s load our data:
##Read in data
#kcal
trainingdata_kcal <- read_excel("ID Training Data copy.xlsx", sheet = 32, range = "A20:I31", col_names = column_names[1:9])
trainingdata_kcal_filtered <- trainingdata_kcal
trainingdata_kcal_filtered[trainingdata_kcal_filtered < 0] <- 0
trainingdata_kcal_filtered_long <- trainingdata_kcal_filtered %>%
pivot_longer(cols = diet_names, names_to = "diet", values_to = "kcal") %>%
subset(select = -c(h2o, spill))
#grams
trainingdata_grams <- read_excel("ID Training Data copy.xlsx", sheet = 32, range = "A3:I14", col_names = column_names[1:9])
trainingdata_grams_filtered <- trainingdata_grams
trainingdata_grams_filtered[trainingdata_grams_filtered < 0] <- 0
trainingdata_grams_filtered_long <- trainingdata_grams_filtered %>%
pivot_longer(cols = diet_names, names_to = "diet", values_to = "grams") %>%
subset(select = -c(h2o, spill))And filter for pure fat and pure sugar:
trainingdata_kcal_filtered_long_pure <- trainingdata_kcal_filtered_long %>%
filter(diet == "zero" | diet == "hundred")
trainingdata_grams_filtered_long_pure <- trainingdata_grams_filtered_long %>%
filter(diet == "zero" | diet == "hundred")And finally we can visualize
trainingdata_kcal_filtered_long_pure %>%
ggplot(mapping = aes(x = rat,
y = kcal,
fill = diet)) +
geom_col(position = "dodge")
trainingdata_grams_filtered_long_pure %>%
ggplot(mapping = aes(x = rat,
y = grams,
fill = diet)) +
geom_col(position = "dodge")
Again, the weight data I think is the best way to go here. Either way, I think a viable strategy here is just to divide rats into fat and sugar preferrers based on the relative amounts they ate here. In that case:
Fat preferrers: ID1, ID12, ID3, ID5, ID9
Sugar preferrers: ID10, ID11, ID2, ID4, ID6, ID7, ID8
Lastly, just to get a sense of each individual rat’s baseline preferences across all diets, I just want to directly plot that as a bar graph:
trainingdata_grams_filtered_long$diet <- factor(trainingdata_grams_filtered_long$diet,
levels = diet_names)
trainingdata_grams_filtered_long$rat <- factor(trainingdata_grams_filtered_long$rat,
levels = rat_names)
trainingdata_grams_filtered_long %>%
ggplot(mapping = aes(x = diet,
y = grams,
fill = diet)) +
geom_col(position = "dodge") +
facet_wrap(~rat)
Based on this, I actually do think a composite measure of baseline preference that accounts of all diets is the best gauge of what the rats actually like. Some of them (like ID12) clearly preferred the highly sugary 80% diet, but didn’t actually eat much of the pure sugar. If we quantified his baseline preferences just from pure sugar, therefore, we’d probably miss out on what his preferences actually were.
Question 8
Part 1
Now let’s try to repeat the analyses from Gosnell et al. In order to repeat their correlation analyses, I need data frames for each non-vehicle diet that contains the following variables:
- Average baseline fat consumption
- Average baseline carbohydrate consumption
- Average baseline for each individual diet
- A column for each diet reflecting the effect of the drug (i.e., the difference between consumption of the that diet for the current drug dose vs. vehicle)
This is going to be hard, but let’s give it a shot.
We already have a dataset containing each rat’s baseline preferences for each diet, so we can start from there by adding composite fat and sugar variables:
baseline_preferences_grams <- trainingdata_grams %>%
subset(select = -c(h2o, spill)) %>%
mutate(
fat = zero + (twenty*0.8) + (forty*0.6) + (sixty*0.4) + (eighty*0.2),
sugar = hundred + (eighty*0.8) + (sixty*0.6) + (forty*0.4) + (twenty*0.2)
)Now we need to make data frames containing intakes for each drug group:
veh_intakes_grams <- data_grams_filtered %>%
filter(condition == "VEH") %>%
subset(select = -c(condition, day)) %>%
rename(veh_zero = zero,
veh_twenty = twenty,
veh_forty = forty,
veh_sixty = sixty,
veh_eighty = eighty,
veh_hundred = hundred,
veh_h2o = h2o,
veh_spill = spill)
hi_intakes_grams <- data_grams_filtered %>%
filter(condition == "HI") %>%
subset(select = -c(condition, day)) %>%
rename(hi_zero = zero,
hi_twenty = twenty,
hi_forty = forty,
hi_sixty = sixty,
hi_eighty = eighty,
hi_hundred = hundred,
hi_h2o = h2o,
hi_spill = spill)
lo_intakes_grams <- data_grams_filtered %>%
filter(condition == "LO") %>%
subset(select = -c(condition, day)) %>%
rename(lo_zero = zero,
lo_twenty = twenty,
lo_forty = forty,
lo_sixty = sixty,
lo_eighty = eighty,
lo_hundred = hundred,
lo_h2o = h2o,
lo_spill = spill)
nal_intakes_grams <- data_grams_filtered %>%
filter(condition == "NAL") %>%
subset(select = -c(condition, day)) %>%
rename(nal_zero = zero,
nal_twenty = twenty,
nal_forty = forty,
nal_sixty = sixty,
nal_eighty = eighty,
nal_hundred = hundred,
nal_h2o = h2o,
nal_spill = spill)And we need to join the vehicle data frame with that of each other drug dose:
hi_effects_grams <- full_join(veh_intakes_grams,
hi_intakes_grams,
by = "rat")
lo_effects_grams <- full_join(veh_intakes_grams,
lo_intakes_grams,
by = "rat")
nal_effects_grams <- full_join(veh_intakes_grams,
nal_intakes_grams,
by = "rat")And now we can calculate drug effects by subtracting the intake of each diet on vehicle with the intake of the same diet on drug:
hi_effects_grams <- hi_effects_grams %>%
mutate(zero_effect = hi_zero - veh_zero,
twenty_effect = hi_twenty - veh_twenty,
forty_effect = hi_forty - veh_forty,
sixty_effect = hi_sixty - veh_sixty,
eighty_effect = hi_eighty - veh_eighty,
hundred_effect = hi_hundred - veh_hundred,
h2o_effect = hi_h2o - veh_h2o,
spill_effect = hi_spill - veh_spill)
lo_effects_grams <- lo_effects_grams %>%
mutate(zero_effect = lo_zero - veh_zero,
twenty_effect = lo_twenty - veh_twenty,
forty_effect = lo_forty - veh_forty,
sixty_effect = lo_sixty - veh_sixty,
eighty_effect = lo_eighty - veh_eighty,
hundred_effect = lo_hundred - veh_hundred,
h2o_effect = lo_h2o - veh_h2o,
spill_effect = lo_spill - veh_spill)
nal_effects_grams <- nal_effects_grams %>%
mutate(zero_effect = nal_zero - veh_zero,
twenty_effect = nal_twenty - veh_twenty,
forty_effect = nal_forty - veh_forty,
sixty_effect = nal_sixty - veh_sixty,
eighty_effect = nal_eighty - veh_eighty,
hundred_effect = nal_hundred - veh_hundred,
h2o_effect = nal_h2o - veh_h2o,
spill_effect = nal_spill - veh_spill)I have a feeling it’ll also make our lives easier if we also join the preference data to each og the effect data frames too:
hi_effects_grams_full <- hi_effects_grams %>%
left_join(baseline_preferences_grams,
by = "rat") %>%
rename(baseline_zero = zero,
baseline_twenty = twenty,
baseline_forty = forty,
baseline_sixty = sixty,
baseline_eighty = eighty,
baseline_hundred = hundred,
baseline_fat = fat,
baseline_sugar = sugar)
lo_effects_grams_full <- lo_effects_grams %>%
left_join(baseline_preferences_grams,
by = "rat") %>%
rename(baseline_zero = zero,
baseline_twenty = twenty,
baseline_forty = forty,
baseline_sixty = sixty,
baseline_eighty = eighty,
baseline_hundred = hundred,
baseline_fat = fat,
baseline_sugar = sugar)
nal_effects_grams_full <- nal_effects_grams %>%
left_join(baseline_preferences_grams,
by = "rat") %>%
rename(baseline_zero = zero,
baseline_twenty = twenty,
baseline_forty = forty,
baseline_sixty = sixty,
baseline_eighty = eighty,
baseline_hundred = hundred,
baseline_fat = fat,
baseline_sugar = sugar)And now we can actually try to run the correlations! Let’s start with the high condition:
hi.rcorr <- rcorr(as.matrix(hi_effects_grams_full[18:33]))
#hi.rcorr
corrplot(hi.rcorr$r)
Now low:
lo.rcorr <- rcorr(as.matrix(lo_effects_grams_full[18:33]))
#lo.rcorr
corrplot(lo.rcorr$r)
And naltrexone:
nal.rcorr <- rcorr(as.matrix(nal_effects_grams_full[18:33]))
#nal.rcorr
corrplot(nal.rcorr$r)
In the interest of not running a billion significance tests for these correlations, I’m going to draw a subset of the ones I’m most interested in and type them into excel. The file is called “Gosnell results”.
The long and short of it is that this analysis doesn’t suggest that baseline preferences have anything to do with the effects of DAMGO However, I’m going to try to get at this question one more way before I give it up.
Let’s imagine you’re a rat in this study, and you’ve been eating these diets for a while. I think it’s pretty unlikely that the stored representation you have of the diets resembles a thought like “I like fat X amount and sugar Y amount, and will therefore eat from all of the diets in weighted proportion to these preferences.” Instead, I think a much more straightforward way of thinking about it would be something closer to “The sixty percent diet is my favorite, so I’m going to beeline to that one.” So, if DAMGO is really working to exacerbate baseline preferences, what you would expect would be an amplification of that effect. Therefore, there should be a positive correlation between intake of each rat’s favorite diet and the effect of morphine on that diet. Similarly, the correlation should be negative for each rat’s least favorite diet.
Part 2
We’ve already calculated and plotted each rat’s baseline preferences:
trainingdata_grams_filtered_long %>%
ggplot(mapping = aes(x = diet,
y = grams,
fill = diet)) +
geom_col(position = "dodge") +
facet_wrap(~rat,
ncol = 3)
And it’s at this point that I realize: almost everyone’s favorite diet
was the 80%. The DAMGO effect being dependent on baseline preferences is
only cool if the rats have different baseline preferences. If
their preferences are all the same, then you would expect exactly what
we’ve seen – higher doses of DAMGO causing increased consumption of that
diet.
This may just be luck of the draw that we got animals who all happened to be 80% preferrers. Alternatively, it may be the case that all rats would prefer 80%. Our dataset doesn’t let us distinguish, but it does reinforce that DAMGO increases consumption of preferred diets.
And just like that, I think we’re done with experiment 2 analysis!
Sike, more plots I want
I want to see each animal’s preferences across testing days too
data_long_foodOnly_factor %>%
ggplot(mapping = aes(x = diet,
y = kcal,
fill = diet)) +
geom_col(position = "dodge") +
facet_wrap(~rat)
data_grams_long_foodOnly_factor %>%
ggplot(mapping = aes(x = diet,
y = grams,
fill = diet)) +
geom_col(position = "dodge") +
facet_wrap(~rat)
I also need to run pairwise comparisons for the one-way ANOVA specifically of the filtered naltrexone data
#Significance
data_grams_long_foodOnly_factor_nal %>%
pairwise_t_test(grams~diet, paired = TRUE,
p.adjust.method = "bonferroni",
estimate = TRUE)## # A tibble: 15 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 grams zero twenty 11 11 0.125 10 0.903 1 ns
## 2 grams zero forty 11 11 0.994 10 0.344 1 ns
## 3 grams zero sixty 11 11 -1.12 10 0.29 1 ns
## 4 grams zero eighty 11 11 -2.64 10 0.025 0.374 ns
## 5 grams zero hundred 11 11 1.33 10 0.213 1 ns
## 6 grams twenty forty 11 11 0.947 10 0.366 1 ns
## 7 grams twenty sixty 11 11 -1.19 10 0.26 1 ns
## 8 grams twenty eighty 11 11 -2.89 10 0.016 0.243 ns
## 9 grams twenty hundred 11 11 1.07 10 0.31 1 ns
## 10 grams forty sixty 11 11 -1.64 10 0.132 1 ns
## 11 grams forty eighty 11 11 -4.11 10 0.002 0.032 *
## 12 grams forty hundred 11 11 1.03 10 0.329 1 ns
## 13 grams sixty eighty 11 11 -0.666 10 0.521 1 ns
## 14 grams sixty hundred 11 11 1.76 10 0.109 1 ns
## 15 grams eighty hundred 11 11 4.42 10 0.001 0.019 *I need to reformat some data for easier use in SPSS
data_kcal_wide <- data_long_foodOnly %>%
subset(select = -day) %>%
pivot_wider(names_from = c(condition, diet), values_from = kcal)
#write_xlsx(data_kcal_wide, "/Users/lawilson1999/Library/Mobile Documents/com~apple~CloudDocs/Psychology MA/Second Year (2022-2023)/Thesis/Thesis analyses (ID)\\data_kcal_wide.xlsx")
data_grams_wide <- data_grams_long_foodOnly %>%
subset(select = -day) %>%
pivot_wider(names_from = c(condition, diet), values_from = grams)
#write_xlsx(data_kcal_wide, "/Users/lawilson1999/Library/Mobile Documents/com~apple~CloudDocs/Psychology MA/Second Year (2022-2023)/Thesis/Thesis analyses (ID)\\data_grams_wide.xlsx")I want to run additional tests comparing the effect of drug dose across diet groups for kcals:
###Filter out diet groups
#kcals
# 0% dataset
data_long_foodOnly_factor_zero <- data_long_foodOnly_factor %>%
filter(diet == "zero")
#20% dataset
data_long_foodOnly_factor_twenty <- data_long_foodOnly_factor %>%
filter(diet == "twenty")
#40%
data_long_foodOnly_factor_forty <- data_long_foodOnly_factor %>%
filter(diet == "forty")
#60%
data_long_foodOnly_factor_sixty <- data_long_foodOnly_factor %>%
filter(diet == "sixty")
#80%
data_long_foodOnly_factor_eighty <- data_long_foodOnly_factor %>%
filter(diet == "eighty")
#100%
data_long_foodOnly_factor_hundred <- data_long_foodOnly_factor %>%
filter(diet == "hundred")
##ANOVAs
#0%
zero.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_zero,
dv = .(kcal),
wid = .(rat),
within = .(condition),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.zero.kcal.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 condition 3 30 1.436595 0.2516532 0.0223938#20%
twenty.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_twenty,
dv = .(kcal),
wid = .(rat),
within = .(condition),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.twenty.kcal.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 condition 3 30 0.9382934 0.434407 0.01115226#40%
forty.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_forty,
dv = .(kcal),
wid = .(rat),
within = .(condition),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.forty.kcal.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 condition 3 30 1.086801 0.369725 0.06477213#60%
sixty.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_sixty,
dv = .(kcal),
wid = .(rat),
within = .(condition),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.sixty.kcal.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 condition 3 30 2.580396 0.07197252 0.08632641#80%
eighty.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_eighty,
dv = .(kcal),
wid = .(rat),
within = .(condition),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.eighty.kcal.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 condition 3 30 9.205299 0.0001795536 * 0.3107644pairwiseEighty <- data_long_foodOnly_factor_eighty %>%
pairwise_t_test(kcal~condition, paired = TRUE,
p.adjust.method = "bonferroni",
estimate = TRUE)
pairwiseEighty## # A tibble: 6 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 kcal HI LO 11 11 2.08 10 0.065 0.388 ns
## 2 kcal HI VEH 11 11 2.28 10 0.046 0.277 ns
## 3 kcal HI NAL 11 11 3.14 10 0.01 0.062 ns
## 4 kcal LO VEH 11 11 0.554 10 0.592 1 ns
## 5 kcal LO NAL 11 11 1.52 10 0.159 0.954 ns
## 6 kcal VEH NAL 11 11 1.99 10 0.075 0.448 ns#100%
hundred.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_hundred,
dv = .(kcal),
wid = .(rat),
within = .(condition),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.hundred.kcal.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 condition 3 30 1.683584 0.1915543 0.1022568And for grams
###Filter out diet groups
#kcals
# 0% dataset
data_grams_long_foodOnly_factor_zero <- data_grams_long_foodOnly_factor %>%
filter(diet == "zero")
#20% dataset
data_grams_long_foodOnly_factor_twenty <- data_grams_long_foodOnly_factor %>%
filter(diet == "twenty")
#40%
data_grams_long_foodOnly_factor_forty <- data_grams_long_foodOnly_factor %>%
filter(diet == "forty")
#60%
data_grams_long_foodOnly_factor_sixty <- data_grams_long_foodOnly_factor %>%
filter(diet == "sixty")
#80%
data_grams_long_foodOnly_factor_eighty <- data_grams_long_foodOnly_factor %>%
filter(diet == "eighty")
#100%
data_grams_long_foodOnly_factor_hundred <- data_grams_long_foodOnly_factor %>%
filter(diet == "hundred")
##ANOVAs
#0%
zero.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_zero,
dv = .(grams),
wid = .(rat),
within = .(condition),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.zero.grams.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 condition 3 30 1.436595 0.2516532 0.0223938#20%
twenty.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_twenty,
dv = .(grams),
wid = .(rat),
within = .(condition),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.twenty.grams.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 condition 3 30 0.9382934 0.434407 0.01115226#40%
forty.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_forty,
dv = .(grams),
wid = .(rat),
within = .(condition),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.forty.grams.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 condition 3 30 1.086801 0.369725 0.06477213#60%
sixty.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_sixty,
dv = .(grams),
wid = .(rat),
within = .(condition),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.sixty.grams.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 condition 3 30 2.580396 0.07197252 0.08632641#80%
eighty.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_eighty,
dv = .(grams),
wid = .(rat),
within = .(condition),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.eighty.grams.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 condition 3 30 9.205299 0.0001795536 * 0.3107644#100%
hundred.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_hundred,
dv = .(grams),
wid = .(rat),
within = .(condition),
return_aov = TRUE,
type = 2)## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.hundred.grams.aov$ANOVA## Effect DFn DFd F p p<.05 ges
## 2 condition 3 30 1.683584 0.1915543 0.1022568Scatters showing correlation between preference and intake:
hi_effects_grams_full %>%
ggplot(mapping = aes(x = baseline_eighty,
y = eighty_effect)) +
geom_point() +
geom_smooth(method = "lm") +
labs(title = "High Dose (0.25µg/0.5µL/side)",
x = "Baseline 80% Sugar Intake (grams)",
y = "Effect of DAMGO on 80% Sugar Intake (grams)") +
theme_classic(base_size = 16)## `geom_smooth()` using formula = 'y ~ x'
lo_effects_grams_full %>%
ggplot(mapping = aes(x = baseline_eighty,
y = eighty_effect)) +
geom_point() +
geom_smooth(method = "lm", se = 0.95) +
labs(title = "Low Dose (0.025µg/0.5µL/side)",
x = "Baseline 80% Sugar Intake (grams)",
y = "Effect of DAMGO on 80% Sugar Intake (grams)") +
theme_classic(base_size = 16)## `geom_smooth()` using formula = 'y ~ x'
Let’s try to get these scatterplots on the same plot:
hi_effects_eighty <- hi_effects_grams_full %>%
subset(select = c(rat, baseline_eighty, eighty_effect)) %>%
mutate(dose = "HI")
lo_effects_eighty <- lo_effects_grams_full %>%
subset(select = c(rat, baseline_eighty, eighty_effect)) %>%
mutate(dose = "LO")
damgo_effects <- rbind(hi_effects_eighty, lo_effects_eighty)
damgo_effects %>%
ggplot(mapping = aes(x = baseline_eighty,
y = eighty_effect,
color = dose)) +
geom_point() +
geom_smooth(method = "lm") +
labs(
x = "Baseline 80% Sugar Intake (grams)",
y = "Effect of DAMGO on 80% Sugar Intake (grams)") +
theme_classic(base_size = 13) +
scale_color_discrete(name = "DAMGO Dose", labels = c("High (0.25 µg/0.5µL/side)", "Low (0.025µg/0.5µL/side"))## `geom_smooth()` using formula = 'y ~ x'